
%----------------------------------------------------------------
% Options
%----------------------------------------------------------------

% @#define alternative_taylor_rule=0 


%----------------------------------------------------------------
% Define variables and parameters
%----------------------------------------------------------------

var pi p_1 p_2 p_3 p_4 p_5
    s mctilde y c n tfp
    markup_1 markup_2 markup_3 markup_4 markup_5
    i a nu i_ann pi_ann y_nat y_nat_mis;

varexo  eps_i
        eps_a
        ;   

parameters  betta siggma varphi eta theta_1 theta_2 theta_3 theta_4 theta_5 phi_pi phi_y  
            rho_a rho_i sig_eps_i sig_eps_a  logpibar  
            logpbar_1 logpbar_2 logpbar_3 logpbar_4 logpbar_5 logsbar
            logybar logybar_sss logmctildebar 
            logmarkupbar_1 logmarkupbar_2 logmarkupbar_3 logmarkupbar_4 logmarkupbar_5 
            logabar lognbar logtfpbar logibar 
            posmp;


%----------------------------------------------------------------
% Parametrization
%----------------------------------------------------------------

load paramfile; 
set_param_value('betta',betta)   
set_param_value('siggma',siggma)  
set_param_value('varphi',varphi)  
set_param_value('eta',eta)     
set_param_value('theta_1',theta_1)  
set_param_value('theta_2',theta_2)  
set_param_value('theta_3',theta_3)  
set_param_value('theta_4',theta_4)  
set_param_value('theta_5',theta_5)  
set_param_value('phi_pi',phi_pi)  
set_param_value('phi_y',phi_y)   

set_param_value('rho_a',rho_a)   
set_param_value('rho_i',rho_i)   

set_param_value('sig_eps_i',sig_eps_i)  
set_param_value('sig_eps_a',sig_eps_a)  

set_param_value('logpibar',log(pibar))   
set_param_value('logpbar_1',log(pbar_1))  
set_param_value('logpbar_2',log(pbar_2)) 
set_param_value('logpbar_3',log(pbar_3)) 
set_param_value('logpbar_4',log(pbar_4)) 
set_param_value('logpbar_5',log(pbar_5)) 
set_param_value('logybar',log(ybar)) 
set_param_value('logybar_sss',log(ybar_sss)) 
set_param_value('logmctildebar',log(mctildebar)) 
set_param_value('logsbar',log(sbar)) 
set_param_value('logmarkupbar_1',log(markupbar_1)) 
set_param_value('logmarkupbar_2',log(markupbar_2)) 
set_param_value('logmarkupbar_3',log(markupbar_3)) 
set_param_value('logmarkupbar_4',log(markupbar_4))
set_param_value('logmarkupbar_5',log(markupbar_5))  
set_param_value('logabar',log(abar)) 
set_param_value('lognbar',log(nbar)) 
set_param_value('logtfpbar',log(tfpbar))  
set_param_value('logibar',log(ibar))   

set_param_value('posmp',posmp)   


%----------------------------------------------------------------
% Equilibrium conditions
%----------------------------------------------------------------

model; 
    ((1-eta) * exp(p_1) + eta * exp(mctilde)) * exp(p_1) ^ (-eta) * exp(y) = theta_1 * (exp(p_1) / exp(p_1(-1)) * exp(pi) - 1) * exp(p_1) / exp(p_1(-1)) * exp(pi) - theta_1 * (exp(p_1(+1)) / exp(p_1) * exp(pi(+1)) - 1) * exp(p_1(+1)) / exp(p_1) * exp(pi(+1));
    ((1-eta) * exp(p_2) + eta * exp(mctilde)) * exp(p_2) ^ (-eta) * exp(y) = theta_2 * (exp(p_2) / exp(p_2(-1)) * exp(pi) - 1) * exp(p_2) / exp(p_2(-1)) * exp(pi) - theta_2 * (exp(p_2(+1)) / exp(p_2) * exp(pi(+1)) - 1) * exp(p_2(+1)) / exp(p_2) * exp(pi(+1));
    ((1-eta) * exp(p_3) + eta * exp(mctilde)) * exp(p_3) ^ (-eta) * exp(y) = theta_3 * (exp(p_3) / exp(p_3(-1)) * exp(pi) - 1) * exp(p_3) / exp(p_3(-1)) * exp(pi) - theta_3 * (exp(p_3(+1)) / exp(p_3) * exp(pi(+1)) - 1) * exp(p_3(+1)) / exp(p_3) * exp(pi(+1));
    ((1-eta) * exp(p_4) + eta * exp(mctilde)) * exp(p_4) ^ (-eta) * exp(y) = theta_4 * (exp(p_4) / exp(p_4(-1)) * exp(pi) - 1) * exp(p_4) / exp(p_4(-1)) * exp(pi) - theta_4 * (exp(p_4(+1)) / exp(p_4) * exp(pi(+1)) - 1) * exp(p_4(+1)) / exp(p_4) * exp(pi(+1));
    ((1-eta) * exp(p_5) + eta * exp(mctilde)) * exp(p_5) ^ (-eta) * exp(y) = theta_5 * (exp(p_5) / exp(p_5(-1)) * exp(pi) - 1) * exp(p_5) / exp(p_5(-1)) * exp(pi) - theta_5 * (exp(p_5(+1)) / exp(p_5) * exp(pi(+1)) - 1) * exp(p_5(+1)) / exp(p_5) * exp(pi(+1));

    1 = 1/5 * exp(p_1) ^ (1-eta) + 1/5 * exp(p_2) ^ (1-eta) + 1/5 * exp(p_3) ^ (1-eta) + 1/5 * exp(p_4) ^ (1-eta) + 1/5 * exp(p_5) ^ (1-eta);

    exp(s) = 1/5 * exp(p_1) ^ (-eta) + 1/5 * exp(p_2) ^ (-eta) + 1/5 * exp(p_3) ^ (-eta) + 1/5 * exp(p_4) ^ (-eta) + 1/5 * exp(p_5) ^ (-eta);

    exp(mctilde) =  exp(s) ^ varphi * exp(c) ^ siggma * exp(y) ^ varphi * exp(a) ^ (-(1+varphi));

    exp(c) ^ (-siggma) = betta * exp(c(+1)) ^ (-siggma) * exp(i) / exp(pi(+1));

    exp(y) = exp(c) + theta_1 * (exp(p_1) / exp(p_1(-1)) * exp(pi) - 1)^2 + theta_2 * (exp(p_2) / exp(p_2(-1)) * exp(pi) - 1)^2 + theta_3 * (exp(p_3) / exp(p_3(-1)) * exp(pi) - 1)^2 + theta_4 * (exp(p_4) / exp(p_4(-1)) * exp(pi) - 1)^2 + theta_5 * (exp(p_5) / exp(p_5(-1)) * exp(pi) - 1)^2;

    exp(n) = exp(y) / exp(a) * exp(s);

    exp(tfp) = exp(y) / exp(n);

    exp(markup_1) = exp(p_1) / exp(mctilde);
    exp(markup_2) = exp(p_2) / exp(mctilde);
    exp(markup_3) = exp(p_3) / exp(mctilde);
    exp(markup_4) = exp(p_4) / exp(mctilde);
    exp(markup_5) = exp(p_5) / exp(mctilde);
 

    @#if alternative_taylor_rule==0
        exp(i) = rho_i * exp(i(-1)) + (1 - rho_i) * (exp(logibar) * (exp(pi) / exp(logpibar)) ^ phi_pi * (exp(y) / exp(y_nat)) ^ phi_y) + nu;
    @#endif
    @#if alternative_taylor_rule==1
        exp(i) = rho_i * exp(i(-1)) + (1 - rho_i) * (exp(logibar) * (exp(pi) / exp(logpibar)) ^ phi_pi * (exp(y) / exp(y_nat_mis)) ^ phi_y) + nu;
    @#endif

    nu = posmp * sig_eps_i * eps_i;

    a  = rho_a * a(-1) + sig_eps_a * eps_a;

    i_ann = 4*i;
    pi_ann = 4*pi;

    exp(y_nat)^(-(siggma+varphi))*exp(a)^(1+varphi) = eta/(eta-1);
    exp(y_nat_mis)^(-(siggma+varphi))*exp(tfp)^(1+varphi) = eta/(eta-1);
end; 

 
%----------------------------------------------------------------
%  steady states
%---------------------------------------------------------------- 

initval;
p_1 = logpbar_1;
p_2 = logpbar_2;
p_3 = logpbar_3;
p_4 = logpbar_4;
p_5 = logpbar_5;
pi = logpibar;
y = logybar;
n = lognbar;
mctilde = logmctildebar;
markup_1 = logmarkupbar_1;
markup_2 = logmarkupbar_2;
markup_3 = logmarkupbar_3;
markup_4 = logmarkupbar_4;
markup_5 = logmarkupbar_5;
c = logybar;
a = logabar;
i = logibar;
s = logsbar;
eps_i = 0;
tfp = logtfpbar;
eps_a = 0;
eps_i = 0;
nu = 0;
y_nat = logybar;
% i_ann = logi_annbar;
% pi_ann = pi_annbar;
end; 

steady;
check;


%----------------------------------------------------------------
% Shocks
%----------------------------------------------------------------

shocks;
    var eps_i = 1;   
    var eps_a = 1;   
end; 



%----------------------------------------------------------------
% Solve
%----------------------------------------------------------------

if order==1
    stoch_simul(order=1, irf=17, irf_shocks=(eps_i, eps_a), nomoments, nodecomposition, nocorr, nofunctions, noprint) pi tfp y i_ann n mctilde s p_1 p_5 markup_1 markup_5; 
elseif order==2
    stoch_simul(order=2, pruning, irf=17, irf_shocks=(eps_i), nomoments, nodecomposition, nocorr, nofunctions, noprint) pi tfp y i_ann n mctilde s p_1 p_5 markup_1 markup_2 markup_3 markup_4 markup_5; 
elseif order==3 
    stoch_simul(order=3, irf=17, irf_shocks=(eps_i, eps_a), nomoments, nodecomposition, nocorr, nofunctions, noprint) pi tfp y i_ann markup_1 n mctilde; 
end
